é a dispersão, desvio ou avermelhamento (ou todos esses efeitos) causados na luz emitida por uma fonte devido ao meio interstelar.
A extinção (ou reddening) é simbolizado pela letra $A_\lambda$ com um subscrito que identifica a banda da correção.
Existem alguns mapas de extinção, sendo Zaritsky, (2004) um dos mais citados. Porém, vou utilizar um mapa criado a partir dos dados de RRLyr da LMC contidos no OGLE-III (Pejcha, 2009) e disponíveis em http://www.astro.princeton.edu/~pejcha/lmc_extmap/.
O arquivo lmcextmap.dat possui os dados do mapa de extinção na seguinte ordem:
MinX, MaxX, MinY, MaxY, E(V-I), E_rms, Delta E(V-I), numero de estrelas
em que as coordenadas $x$ e $y$ são obtidas pela conversão, $$ x = (\alpha -80.35)\cdot \cos(\delta) $$ $$ y = \delta +69.68 $$ onde $\alpha$ é ascenção reta e $\delta$ é a declinação, ambas dadas em graus.
Para obter a extinção, devo converter os dados do OGLE para as coordenadas $x$ e $y$ e obter E(V-I) da tabela lmcextmap.dat. A extinção na banda I é obtida pelo calculo, $$A_I = R_{VI} E(V-I)$$ em que $R_{VI} = 1.10$ como utilizado no paper.
Porem, no catálogo OGLE, $\alpha$ é dado em horas. Primeiramente vou converter $\alpha$ em graus.
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
#Ler a tabela de dados com RA e Dec de cada RRLyr AB
# a estrela OGLE-LMC-RRLYR-15485 (posição 11013 no arquivo) não possui DatabaseNumber, estava dando problema para ler os dados. Coloquei manualmente o valor 0000
from astropy.io import ascii
tbl = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/identification_onlyRRAB.dat', guess=False, Reader=ascii.Tab, header_start = 0, data_start = 1)
tbl
Out[1]:
In [2]:
#Converter RA de horas para graus
import astropy.coordinates as coord
import astropy.units as u
ra = coord.Angle(tbl['RA'], unit = u.hour)
ra.degree
Out[2]:
In [3]:
#Converter Dec em graus
dec = coord.Angle(tbl['Dec'], unit = u.degree)
dec.degree
Out[3]:
Agora, tendo os valores de RA e Dec em graus, vou converter para o sistema X e Y utilizado no paper
In [4]:
#Converter RA e Dec em X e Y de acordo com o paper do Pejcha (2009)
import numpy as np
x = (ra.degree - 80.35) * np.cos(dec.degree)
y = dec.degree + 69.68
ra_dec = np.around(np.column_stack((x,y)), decimals=1) #arrendondamento para o primeiro decimal apos a virgula
print(ra_dec)
#np.savetxt('ra_dec.txt',ra_dec, fmt='%s')
Ler os dados do extinction map em forma de tabela
In [5]:
#ext_map = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/extinction_map/lmcextmap.dat',
# guess=False, Reader=ascii.Basic)
ext_map = np.array(np.loadtxt('/home/glauffer/Dropbox/FURG/final_project/data/extinction_map/lmcextmap.dat', skiprows=29))
ext_map
Out[5]:
In [6]:
#Criar array 2d com os dados min e max para X e Y
#x_map = np.column_stack((ext_map['minX'], ext_map['maxX']))
#y_map = np.column_stack((ext_map['minY'], ext_map['maxY']))
x_map = np.column_stack((ext_map.T[0], ext_map.T[1]))
y_map = np.column_stack((ext_map.T[2], ext_map.T[3]))
ext = np.ones(len(ra_dec))
for i in range(len(ra_dec)):
if ra_dec[i,0]>(-4.7): # or ():
ext[i] = 0
if ra_dec[i,0]>4.9:
ext[i] = 0
if ra_dec[i,1]>(-3.4): #or
ext[i] = 0
if ra_dec[i,1]>3.2:
ext[i] = 0
else:
try:
ext[i] = ext_map[np.where((x_map.T[0] == ra_dec[i,0]) & (y_map.T[0] == ra_dec[i,1]))][0][4]
except IndexError:
ext[i] = 0
ext
Out[6]:
In [7]:
#np.savetxt('teste_ext.txt', ext, fmt='%s')
#Ler os dados de magnitude das RRLYR AB
#Nome MagI MagV OGLE_period RA DEC CE_period...
magI = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/ID_mag_Period_RA_DE.dat')
magI
Out[7]:
In [8]:
len(magI['col2'])
Out[8]:
In [9]:
#corrigir a extincao
magI_correct = magI['col2'] - 1.10*ext
In [10]:
len(magI_correct)
Out[10]:
In [11]:
#relacao pl -> magI - ext = a*log(P) + B + u
#em que u eh o modulo da distancia, ou seja, quanto se distancia do eixo x
#plot da dispersao
plt.plot(np.log(magI['col7']), magI_correct, 'bo', ms=0.5)
plt.gca().invert_yaxis()
In [12]:
#Least squares regression for PL relation
A = np.column_stack((np.log(magI['col7']), np.ones_like(magI['col7'])))
x, res, rank, s = np.linalg.lstsq(A,magI_correct)
print('m - A= ', x[0], '*log P + ', x[1])
In [13]:
plt.plot(np.log(magI['col7']), magI_correct, 'bo', ms=0.4)
plt.plot(np.log(magI['col7']), x[0]*np.log(magI['col7']) + x[1], 'r-', label='m - A= %f*log P +%f'%(x[0],x[1]))
plt.gca().invert_yaxis()
plt.legend()
plt.xlim(-2.0,0)
Out[13]:
In [14]:
#u = magI_correct - x[0]*np.log(P) - x[1]
u = np.zeros(len(magI_correct))
u = magI_correct - x[0]*np.log(magI['col7'])-x[1]
In [15]:
u
Out[15]:
In [16]:
#passando para distancia -> r = 10**(0.2*(u+5))
r = np.zeros(len(u))
r = 10**(0.2*u+0.2*5)
In [17]:
r
Out[17]:
In [18]:
#sendo a distancia media 50 kpc -> r_lmc = r+50
r_lmc = np.zeros(len(r))
r_lmc = r + 50e3
In [19]:
#r_lmc
#plot das coordenadas RA e DEC
plt.plot(ra.degree, dec.degree, 'bo', ms=1.)
plt.plot(np.mean(ra.degree), np.mean(dec.degree), 'go')
print(np.mean(ra.hour), np.mean(dec.degree))
#plt.plot(ra.radian, dec.radian, 'bo', ms=1.)
#plt.plot(np.mean(ra.radian), np.mean(dec.radian), 'go')
In [20]:
#theta = dec.radian
theta = ra.radian
theta
Out[20]:
In [21]:
ax = plt.subplot(111, polar=True)
#l=17000
#plt.gca(projection='polar')
#plt.plot(2*np.pi*theta[0:l], r_lmc[0:l], 'ro')#dec.degree, 'ro')
#plt.plot(np.sin(theta), r_lmc, 'ro', ms=0.5)
plt.plot(theta, r_lmc, 'bo', ms=1.0)
ax.set_rlim(0, 55e3)
#ax.set_ylim(30,70)
#ax.set_yticks([35, 50, 60, 70]) #['35 kpc', '50 kpc', '60 kpc', '70 kpc']
#plt.thetagrids([20,40,60], labels=None, fmt='%d', frac = 1.1)
#plt.scatter(theta, r, 'ro')
#plt.show()
#plt.gca(projection='polar')
Out[21]:
In [22]:
#Calculando as coordenadas x, y, z
# x = -r_lmc * np.sin(ra - ra_0) * np.cos(dec)
# y = r_lmc * np.sin(dec)*np.cos(dec_0) − r_lmc * np.sin(dec_0) * np.cos(ra - ra_0) * np.cos(dec)
# z = r_lmc_0 - r_lmc * np.sin(dec) * np.sin(dec_0) - r_lmc * np.cos(dec_0)*np.cos(ra) - \ra_0 * np.cos(dec)
x, y, z = np.zeros(len(r_lmc)), np.zeros(len(r_lmc)), np.zeros(len(r_lmc))
r_lmc_0 = 50e3
ra_0 = np.mean(ra.radian)
dec_0 = np.mean(dec.radian)
x = -r_lmc * np.sin(ra.radian - ra_0) * np.cos(dec.radian)
y = r_lmc * np.sin(dec.radian)*np.cos(dec_0) - r_lmc * np.sin(dec_0) * np.cos(ra.radian - ra_0) * np.cos(dec.radian)
z = r_lmc_0 - r_lmc * np.sin(dec.radian) * np.sin(dec_0) - r_lmc * np.cos(dec_0)*np.cos(ra.radian) - ra_0 * np.cos(dec.radian)
In [24]:
from mpl_toolkits.mplot3d import Axes3D
plt.clf()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.plot(x, y, z, 'bo' ,ms=0.2)
#fig.show()
Out[24]:
In [23]: